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OBJECTIVES 

The  primary  objectives  of  this  research  were  to  develop  reduced-order  modeling  method¬ 
ologies  for  control  and  optimization  applications.  This  was  achieved  through  three  specific 
research  goals.  First  of  all,  we  focused  on  improving  the  accuracy  of  reduced-order  model 
simulations  by  including  closure  terms  (that  can  be  implemented  in  practical  applications). 
Secondly,  we  used  sensitivity  analysis  to  improve  the  accuracy  of  models  over  ranges  of 
parameters.  Finally,  we  developed  a  better  methodology  for  using  reduced-order  modeling 
methods  for  feedback  control  calculations.  We  made  significant  advances  in  each  of  these 
three  areas. 

ACCOMPLISHMENTS 


Overview 

Our  research  program  focused  on  reduced-order  models  for  large-scale  systems  with 
an  emphasis  on  fluids  and  application  to  control  and  optimization  problems.  In  this  sum¬ 
mary,  we  will  focus  on  three  results  that  correspond  to  the  three  main  tasks  outlined  in  our 
original  proposal.  In  the  first  section,  we  provide  an  overview  of  our  progress  introduc¬ 
ing  closure  models  for  POD  /  Galerkin  dynamical  systems  based  on  the  advances  in  the 
turbulence  modeling  community.  Results  are  presented  for  Smagorinsky  closure  models 
in  two-dimensional  and  three-dimensional  flows  about  circular  cylinders.  This  includes 
those  corresponding  to  a  multi-level  discretization  approach  to  make  this  closure  approach 
feasible.  In  the  second  section,  we  highlight  the  use  of  sensitivity  analysis  to  improve  the 
accuracy  of  reduced-order  models  over  ranges  in  parameters.  Each  of  these  two  approaches 
(closure  and  sensitivity)  provide  significant  improvements  to  reduced-order  models  mak¬ 
ing  them  more  attractive  as  surrogates  in  optimal  design  problems.  In  the  third  section, 
we  present  a  novel  approach  to  solving  flow  control  problems  by  developing  reduced-order 
models  specifically  for  Lyapunov  and  Chandrasekhar  equations  (rather  than  designing  con¬ 
trollers  for  reduced-order  models). 


Closure  for  POD  Models 


Accurate  modeling  of  complex  flows  using  the  proper  orthogonal  decomposition  (POD) 
requires  closure  terms  that  account  for  the  discarded  (truncated)  basis  functions.  Prior  to 
this  research,  closure  models  were  either  based  on  tuning  model  coefficients  to  match  data 
or  based  on  simplified  turbulence  models.  The  latter  is  due  to  the  desire  to  build  models 
that  can  be  evaluated  rapidly.  However,  as  our  preliminary  studies  have  shown,  introduc¬ 
ing  models  based  on  modern  turbulent  flow  practices  leads  to  more  accurate  reduced-order 
models  at  a  significantly  higher  computational  cost.  A  survey  of  our  closure  modeling 
approaches  appear  in  papers  #7,  10,  11,  18,  24,  39,  and  47.  An  enabling  method  for  im¬ 
plementing  our  closure  models  is  two-level  (or  multi-level)  discretization  of  the  nonlinear 
closure  terms.  The  general  approach  is  to  use  coarser  meshes  to  discretize  the  nonlinear 
closure  terms  involving  the  lower  order  POD  modes.  We  highlight  the  main  results  from 
this  study  using  a  two-dimensional  flow  past  a  cylinder  at  Reynolds  number  Re  =  200  and 
a  three-dimensional  flow  past  a  cylinder  at  Reynolds  number  Re  =  1000  (details  can  be 
found  in  paper  #  47). 

We  approximate  the  flow  velocity  v  using  a  centering  trajectory  v  and  a  linear  combi¬ 
nation  of  POD  modes  s'pan{(f)j\r]=l, 

r 

v(x,t)  ~  vr(x,t)  =  v(x)  +  ^  (, bj(x)a,j(t ). 

j= i 

Substitution  of  this  form  into  the  weak  form  of  the  Navier-Stokes  equations  leads  to  a 
low-order  dynamical  system  of  the  form 

a  =  b  +  Aa  +  aT Ba,  Oj(0)  =  (4>j,  v(-,  0)  —  v(-)). 

However,  when  comparing  coefficients  a,  with  projected  quantities  dj  =  v(-,  t) 

it  is  obvious  that  one  needs  to  model  the  influence  of  the  discarded  POD  modes  {4>r+i-  •  •  •} 
in  the  dynamical  system.  This  is  more  obvious  as  the  Reynolds  number  increases.  In  the 
study  below,  we  use  a  standard  Smagorinsky  model  to  account  for  the  dissipative  effect  of 
the  smaller  scale  structures  in  the  flow.  We  are  currently  preparing  a  study  that  compares  a 
variety  of  closure  models  using  the  multi-level  discretization  approach. 

In  Figure  1,  we  present  the  coefficients  di  and  a3  as  well  as  the  corresponding  coeffi¬ 
cients  using  the  standard  POD  /  Galerkin  approach  for  the  two-dimensional  cylinder  flow 
problem  (labeled  aPOD ).  Note  there  is  a  poor  agreement  between  these  two  models  (inte¬ 
grated  over  1000  shedding  cycles). 

We  use  the  Smagorinsky  closure  model  account  for  the  discarded  modes  and  present 
the  evolution  of  the  first  and  third  coefficients  in  Figure  2.  These  coefficients  have  much 
better  agreement  with  a,  however,  introduce  a  substantial  computational  cost.  This  cost 
can  be  minimized  by  employing  a  multi-level  discretization  approach  (using  a  4x  coarser 
mesh  for  the  nonlinear  Smagorinsky  terms).  As  seen  in  Figure  2,  there  is  little  difference 
in  the  solutions  using  this  modified  closure  term,  but  the  computational  cost  is  reduced  by 
a  factor  of  nearly  20. 

This  study  was  repeated  for  a  more  challenging  flow  field  corresponding  to  three  dimen¬ 
sional  flow  past  a  cylinder  at  Re  =  1000  (the  computational  mesh  is  indicated  in  Figure  3). 


Figure  1 :  Comparison  of  a  (top)  and  aPOD  (bottom)  coefficients 


Figure  2:  Coefficients  of  reduced-order  model  with  Smagorinsky  closure  terms:  original 
mesh  (top),  coarse  mesh  (bottom) 
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Figure  3:  Computational  grid  for  three-dimensional  flow  past  a  cylinder 


A  reduced-order  model  for  r  =  6  was  computed  and  found  similar  improvements  with  the 
addition  of  a  Smagorinsky  closure  model.  A  table  indicating  the  computational  speedup 
of  the  two-level  methods  over  the  one-level  method  as  well  as  the  relative  error  in  the  flow 
field  over  one  shedding  cycle  is  given  below. 

coarsening  level  speedup  factor  error 
“T  I  4.46  x  KT2 

2  5.22  4.52  x  1(T2 

4  24.18  4.73  x  KT2 


Sensitivity  Analysis  of  POD  Modes 

This  research  led  to  an  extensive  study  on  the  benefits  of  using  sensitivity  analysis  to 
improve  the  quality  of  reduced  order  models  from  a  single  simulation.  The  results  of  this 
research  is  contained  in  references  #5,  9,  25,  30,  31,  32,  33,  34,  35,  and  42.  We  consider 
a  POD  model  for  a  flow  with  a  geometric  parameter  change.  The  flow  again  consists 
of  two  dimensional  flow  past  a  square  cylinder,  but  the  angle  of  the  square  to  the  fixed 
incoming  flow  and  channel  walls  is  parameter  dependent.  First  of  all,  we  introduce  a  mesh 
warping  function  to  transform  the  domain  at  7  =  22.5°,  a  =  0°  to  a  range  of  parameters 
a  (see  Figure  4).  This  mesh  warping  function  is  shown  for  a  =  —22.5°,  0°,  and  22.5°, 
respectively  in  Figure  5. 

This  mesh  warping  function  allows  us  to  map  POD  bases  to  different  geometric  con¬ 
figurations  as  well  as  define  Lagrangian  flow  sensitivity  analysis,  to  describe  how  the  flow 
changes  with  changes  to  the  geometry.  This  sensitivity  analysis  can  be  used  to  extrapolate 
the  nominal  POD  basis  computed  at  a  =  0°  to  a  range  of  parameters.  We  also  consider 
building  reduced-order  models  that  incorporate  both  the  original  POD  basis  as  well  as  the 
geometric  sensitivity  of  the  basis  to  a  (an  example  of  a  POD  basis  function  and  the  sensi¬ 
tivity  of  the  POD  basis  function  with  respect  to  a  parameter  are  provided  in  Figure  7). 

We  compute  the  relative  error  (vs.  direct  numerical  simulation)  in  reduced  order  models 
of  order  12  over  one  shedding  period  and  plot  the  results  in  Figure  6.  We  note  that  the 
extrapolated  basis  provides  a  significant  improvement  in  the  reduced  order  models  (nearly 
an  order  of  magnitude  in  relative  error)  at  small  parameter  changes  where  the  extrapolation 
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Figure  4:  Flow  past  a  cylinder  at  varying  angle  of  attack,  a 


Figure  5:  Demonstration  of  mesh  warping  function 


is  expected  to  be  valid.  The  error  in  the  reduced-order  model  with  the  extrapolated  basis  is 
comparable  to  that  obtained  by  projecting  the  CFD  solution  onto  the  basis. 


Figure  6:  Relative  errors  in  reduced-order  models 

We  note  that  the  example  with  significant  geometric  deformation  is  challenging  for 
most  reduced-order  modeling  capabilities.  Combining  sensitivity  analysis  (for  varying  vis¬ 
cosity  and  boundary  conditions)  with  POD  had  more  dramatic  improvement  over  larger 
parametric  variation. 
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(a.)  POD  modes 


(b.)  Sensitivity  of  POD  modes  with  respect  to  a. 


Figure  7:  POD  modes  and  their  geometric  sensitivity  to  a 


A  Reduced  Order  Solver  for  Lyapunov  Equations  with  High  Rank  Matrices 


Consider  the  Lyapunov  equation 


AtP  +  PA  +  Q  =  0  (1) 

where  A  and  Q  are  sparse1,  high  rank  matrices  and  A  is  asymptotically  stable.  We  further 
consider  the  special  case  where  Q  is  symmetric,  leading  to  symmetric  solutions  P.  Further¬ 
more,  we  are  only  interested  in  the  action  of  P  on  one  column  vector  b  (or  a  few  column 
vectors  in  B),  in  other  words,  we  seek  Pb  (or  PB).  Among  many  control  applications, 
one  natural  application  arises  from  Kleinman-Newton  iterations  to  solve  large  scale  Riccati 
equations. 

It  is  well  known  that  the  solution  P  to  (1)  can  be  written  as 

p=  [°°  eATtQeAt  dt . 

Jo 

Thus,  the  action  of  P  on  a  vector  b  can  be  approximated  as  the  solution  to  a  two  point 
boundary  value  problem.  To  see  this,  we  introduce  the  weakly  coupled  state/adjoint  system 

x  =  Ax,  x(0)  =  b  (2) 

— Ar  =  ArAT  +  Qx,  A  (t)  =  0.  (3) 

Since  limT^.oc,  x(r)  =  0,  our  approximation  is  seen  in  the  limiting  solution 

Aoo(-)  =  lim  At(-). 

t  — yoo 

Since  A <*,(£)  =  Px(t),  we  arrive  at  our  desired  product  as  Aoo(0)  =  Px(0)  =  Pb. 

A  straightforward  approach  for  approximating  Pb  would  be  to  numerically  integrate  (2) 
until  the  solution  is  essentially  zero  (guaranteed  by  the  asymptotic  stability  assumption  on 
A).  With  the  state  solution  x  available,  the  adjoint  equation  can  be  integrated  backward 
in  time.  Since  storage  would  likely  be  an  issue,  a  checkpointing  strategy,  eg.  [Griewank 
2000],  could  be  employed  if  needed  to  carry  out  this  approach.  However,  the  cost  of  inte¬ 
grating  equations  (2)-(3),  even  exploiting  the  fact  that  only  products  Ax,  ATA  and  Qx  are 
needed,  severely  limit  the  applicability  of  this  approach. 

The  central  strategy  proposed  in  this  research  is  to  split  the  integration  above  into  two 
pieces.  The  solution  over  the  interval  [0,  T]  is  integrated  accurately  until  model  reduction 
can  be  performed  to  accurately  approximate  the  solution  over  [ T ,  oo).  Then  the  system  over 
[T,  oo)  reduces  to  a  low  order  Sylvester  equation  from  which  we  can  obtain  an  approximate 
final  condition  for  A(T).  With  the  final  condition  we  can  integrate  the  system  back  to 
t  =  0.  Our  reduced  order  solver  is  outlined  below.  Note  that  for  simplicity,  we  consider 
approximating  the  product  Pb  for  a  column  vector  b.  Motivation  and  parameter  selection 
for  the  algorithm  can  be  found  in  papers  #8  and  26. 

'Alternately,  our  methodology  applied  to  the  case  where  A  and  Q  can  be  decomposed  into  basic  opera¬ 
tions  involving  sparse  matrices  such  that  the  products  Ax,  AT  A  and  Qx  can  be  computed  cheaply 


Algorithm  1  (Reduced  Order  Lyapunov  Solver)  Given  stable  A  e  IRnxn,  b  e  IRnXl, 
Q  G  IRnXn  (Q  =  QT)  as  well  as  parameters  T  and  r  -C  n. 

1.  Integrate  the  state  equation  (2)  from  x(0)  =  b  over  the  interval  [0,  T], 

2.  Use  the  proper  orthogonal  decomposition,  or  another  reduced  order  basis  method, 
to  generate  a  low  r-dimensional  basis  for  x(i)  on  the  time  interval  [T,  oo).  This  basis 
is  arranged  in  columns  of  the  matrix  V  €  IRnxr. 

3.  Solve  the  Sylvester  equation 

AtP  +  PVtAV  +  QV  =  0. 

Note  that  this  can  be  done  by  performing  Schur  decomposition  to  the  low  order  ma¬ 
trix  VTAV  and  then  solving  r  linear  systems  of  order  n  (Bartels  and  Stewart). 

4.  Integrate  the  adjoint  equation  (3)  from  A (T)  =  PVrx(T)  over  [0,  T]. 

The  result  is  A(0)  ~  Pb. 

This  algorithm  was  tested  on  a  two-dimensional  advection  diffusion  reaction  equation  re¬ 
sulting  from  a  linearized  2D  Burgers  equation  [Camphouse  2004],  In  the  table  below,  we 
present  relative  errors  of  the  solution  for  varying  values  of  integration  time  T  and  reduced- 
order  model  size  r. 


Mesh  size  N  =  1413. 

r 

5 

10 

20 

T  =  2 
T  =  4 

T  =  6 

1.5882e-05 

1.2096e-08 

6.5992e-12 

2.0714e-09 

4.6095e-13 

4.6319e-13 

4.6295e-13 

4.6333e-13 

4.6317e-13 

An  Efficient  Long-Time  Integrator  for  Chandrasekhar  Equations 

The  Chandrasekhar  equations  have  been  posed  as  a  methodology  for  solving  the  infinite 
horizon  control  problem  when  the  system  involves  a  distributed  parameter  system.  In  this 
case,  the  Chandrasekhar  equations  have  the  form 

-K(f)  =  R“1BTL(f)LT(f),  K(0)  =  0  G  IRmxn  (4) 

— L (t)  =  (A  -  BK(f))T  L (t),  L(0)  =  CT  e  KnXp.  (5) 

The  solution  to  the  regulator  problem  (A,  B,  C)  is  then  given  by 

K  =  lim  K (t). 

t — y — oo 

This  approach  replaces  the  need  to  find  the  (dense)  n  x  n  solution  to  the  Riccati  equation 
by  the  integration  of  ( m+p)n  equations  towards  a  steady  state  solution.  While  the  reduced 
storage  costs  of  the  Chandrasekhar  equations  make  some  large  problems  tractable  that  may 


not  be  otherwise  solvable,  the  slow  convergence  towards  a  steady  state  solution  magnify 
the  computational  costs  associated  with  the  integration. 

As  with  the  Lyapunov  equation  above,  we  have  developed  an  algorithm  that  uses  reduced- 
order  modeling  ideas  to  dramatically  reduce  the  integration  time  that  is  required  for  (4)-(5). 
Consider  the  expression  (from  the  derivation  of  the  Chandrasekhar  equations) 

II  =  [°  L(t)LT(t)  dt  +  [ T  L(t)LT(t)  dt .  (6) 

J  T  J — oo 

N - v - ' 

Ilres 

The  first  term  above  is  integrated  using  the  Chandrasekhar  equation,  while  we  approximate 
the  second  integral  using  a  low-dimensional  basis  for  L  (over  (— oo,  T)).  The  expression 
(6)  above  shows  a  clear  connection  between  finding  a  good  basis  for  L  and  a  good  basis  for 

hires- 

Algorithm  2  (Long-Time  Integrator  for  (4)-(5))  Given  A,  B  and  C. 

1.  Integrate  the  Chandrasekhar  equations  until  time  T  <  0. 

2.  Build  a  reduced-basis  V  for  L (£),  t  <  T. 

3.  Solve  a  reduced  Riccati  equation  for  P,  set  Kres  =  R_1BTVPVT. 

4.  Use  Kres  =  R_1BTIIres  to  compute  K  =  K(T)  +  Kres. 

The  accuracy  of  Algorithm  2  (relative  error  Ere \  for  different  integration  times  |T|  and 
different  sizes  of  the  reduced-order  model  (r)  was  performed.  As  we  expect,  Ere \  is  re¬ 
duced  with  longer  integration  time  and  larger  model  dimension.  However,  we  observe 
that  longer  integration  times  lead  to  smaller  model  dimension  requirements  to  get  “maxi¬ 
mum”  accuracy.  This  supports  our  rationale  that  components  of  higher  frequency  modes 
are  ultimately  eliminated  with  integration.  It  is  impressive  that  the  relative  error  of  2.8% 
obtained  at  T  =  —10  can  be  reduced  5  orders  of  magnitude  with  a  3  dimensional  model. 
Comparisons  are  provided  in  references  #27  and  46. 


Continuing  Research 


The  co-PIs,  Jeff  Borggaard  and  Traian  Iliescu  along  with  John  Paul  Roop  at  North  Car¬ 
olina  A&  T  are  currently  extending  this  research  toward  the  rapid  simulation  of  Boussinesq 
equations.  This  includes  application  of  the  above  techniques  with  the  goal  of  improving  the 
design  /  control  of  thermal  fluids  in  buildings  as  well  as  other  applications  in  geophysical 
fluids. 

Natural  extensions,  such  as  Hermite  interpolation  of  POD  bases  in  parameter  space, 
incorporating  other  subgridscale  and  dynamic  LES  models,  and  their  combination  are  cur¬ 
rently  underway.  Continuation  of  the  promising  flow  control  approaches  discovered  in  this 
research  are  also  planned. 
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2011). 

2.  201 1  SIAM  Student  Conference,  Clemson  University  (February  201 1). 

3.  The  30th  Southeastem-Atlantic  Regional  Conference  on  Differential  Equations,  Vir¬ 
ginia  Tech  (October  2010). 

4.  34th  SIAM  Southeastern- Atlantic  Section  Conference,  North  Carolina  State  Univer¬ 
sity,  (March  2010). 

5.  SIAM  Student  Conference  2010,  Virginia  Tech  (February  2010). 

6.  The  Clemson/Pitt/UTK/VT  graduate/post  graduate  SIAM  Student  conference,  Vir¬ 
ginia  Tech,  Blacksburg,  VA  (February  2009). 

H.-W.  van  Wyk 

1.  35th  SIAM  Southeastern- Atlantic  Section  Conference,  Charlotte,  NC,  (March  201 1). 

2.  2011  SIAM  Student  Conference,  Clemson  University  (February  2011). 

3.  SIAM  Student  Conference  2010,  Virginia  Tech  (February  2010). 

4.  The  Clemson/Pitt/UTK/VT  graduate/post  graduate  SIAM  Student  conference,  Vir¬ 
ginia  Tech,  Blacksburg,  VA  (February  2009). 

Z.  Wang 

1.  2011  SIAM  Computational  Science  and  Engineering,  Reno,  NV  (March  2011). 

2.  201 1  SIAM  Student  Conference,  Clemson  University  (February  201 1). 

3.  The  30th  Southeastern-Atlantic  Regional  Conference  on  Differential  Equations,  Vir¬ 
ginia  Tech  (October  2010). 

4.  Student  Argonne  Summer  Symposium,  Argonne  National  Laboratory  (August  2010). 

5.  2010  SIAM  Annual  Meeting  (AN10),  Pittsburgh,  PA  (July  2010). 

6.  34th  SIAM  Southeastern-Atlantic  Section  Conference,  North  Carolina  State  Univer¬ 
sity,  (March  2010). 

7.  SIAM  Student  Conference  2010,  Virginia  Tech  (February  2010). 

8.  The  First  VT  Symposium  on  Reduced-Order  Modeling  and  System  Identification, 
Virginia  Tech  (February  2010). 


9.  Fall  Fluid  Mechanics  Symposium,  Virginia  Tech  (November  2009). 


10.  SIAM  Conference  on  Computational  Science  &  Engineering  (CSE09),  Miami,  FL 
(March  2009). 

11.  The  Clemson/Pitt/UTK/VT  graduate/post  graduate  SIAM  Student  conference,  Vir¬ 
ginia  Tech,  Blacksburg,  VA  (February  2009). 

12.  Project/NExt/Young  Mathematician’s  Network  Poster  Session  in  AMS  Joint  Mathe¬ 
matics  Meeting,  Washington,  DC  (January  2009). 

Interactions  and  Transitions 

Air  Force  Research  Laboratory,  Wright-Patterson  Air  Force  Base,  OH 

Jeff  Borggaard  spent  the  summer  of  2007  at  AFRL/VACA  working  with  Chris  Camphouse 
and  James  Myatt  on  a  flow  control  problem.  Efficient  POD  software  using  algorithms 
developed  in  this  research  and  applicable  to  practical  engineering  reduced-order  model- 
based  control  algorithms  was  provided  to  the  AFRL  researchers.  Imran  Akhtar  visited 
several  AFRL  researchers  during  a  visit  to  Dayton  in  the  spring  of  2010.  This  included 
several  meetings  with  Phil  Beran  both  in  Dayton  and  when  Dr.  Beran  visited  Virginia 
Tech. 

A  Ph.D.  student  at  the  Interdisciplinary  Center  for  Applied  Mathematics,  Capt.  Kevin 
Pond,  started  working  at  AFIT  in  the  fall  of  2010.  Jeff  Borggaard  visited  AFIT  and  AFRL 
in  the  fall  of  2010,  including  a  meeting  with  Jack  Benek  to  foster  collaborations. 

Honors/Awards 

Jeff  Borggaard  was  awarded  an  ASEE  Summer  Faculty  Fellowship  and  spent  May- July 
2007  at  the  AFRL  Control  Sciences  Center  of  Excellence.  This  occurred  just  prior  to  this 
grant  period. 


